#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/blas_lapack.h>
#include <src/random.h>
#include <src/hf.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   constructs the one-body reduced density matrix from the orbital coefficient matrix
///
/// \param[in]     n        number of orbitals
/// \param[in]     np       number of particles
/// \param[in]     c        the matrix of orbital coefficients,  size n*n
/// \param[out]    x        on exit, the one body reduced density matrix, size n*n
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::build_rdm1(const int n,
                                                  const int np,
                                                  const std::vector<S> & c,
                                                        std::vector<S> & x)
{
  for (int p = 0; p < n; p++)
  for (int q = 0; q < n; q++) {
    x[n*p+q] = pcps::zero<S>();
    for (int i = 0; i < np; i++)
      x[n*p+q] += pcps::conj( c[n*i+p] ) * c[n*i+q];
  }
}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   constructs the alpha fock operator
///
/// \param[in]     n        number of orbitals
/// \param[in]     ra       the one body reduced density matrix for alpha, size n*n
/// \param[in]     rb       the one body reduced density matrix for beta , size n*n
/// \param[in]     oei      the one body integrals, size n*n
/// \param[in]     tei      the two body integrals, which are assumed real, in 1122 order, size n*n*n*n
/// \param[out]    f        on exit, the matrix for the fock operator, size n*n
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::build_fock(const int n,
                                                  const std::vector<S> & ra,
                                                  const std::vector<S> & rb,
                                                  const std::vector<S> & oei,
                                                  const std::vector<S> & tei,
                                                        std::vector<S> & f)
{

  // contribution from the one-electron part of the Hamiltonian
  pcps::xcopy(n*n, &oei.at(0), 1, &f.at(0), 1);

  // coulomb contribution from alpha
  pcps::xgemm('N', 'N', n*n, 1, n*n, pcps::unity<S>(), &tei.at(0), n*n, &ra.at(0), n*n, pcps::unity<S>(), &f.at(0), n*n);

  // coulomb contribution from beta
  pcps::xgemm('N', 'N', n*n, 1, n*n, pcps::unity<S>(), &tei.at(0), n*n, &rb.at(0), n*n, pcps::unity<S>(), &f.at(0), n*n);

  // exchange contribution
  for (int p = 0; p < n; p++)
  for (int q = 0; q < n; q++)
  for (int r = 0; r < n; r++)
  for (int s = 0; s < n; s++)
    f[n*r+q] -= tei[n*n*n*p+n*n*q+n*r+s] * ra[n*p+s];

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   diagonalizes the fock matrix and extracts the new orbital coefficients
///
/// \param[in]     n        number of orbitals
/// \param[in]     np       number of particles
/// \param[in]     f        the matrix for the fock operator, size n*n
/// \param[out]    eval     on exit, the eigenvalues of the the fock matrix, size n
/// \param[out]    c        on exit, the matrix of orbital coefficients, size n*n
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::diagonalize_fock(const int n,
                                                        const int np,
                                                        const std::vector<S> & f,
                                                        std::vector<double> & eval,
                                                        std::vector<S> & c)
{

  // preparation for diagonalization
  int info;
  const int lwork = 5*n;
  std::vector<S> work(lwork);
  std::vector<double> rwork(3*n);

  // move the fock matrix into the orbital coefficient matrix in preparation for diagonalization
  pcps::xcopy(n*n, &f.at(0), 1, &c.at(0), 1);

  // take the complex conjugate of the fock matrix.  This way when FORTRAN transposes it the
  // resulting matrix will be equal to the original fock matrix
  pcps::xconj(n*n, &c.at(0), 1);

  // diagonalize the fock matrix, after which the eigenvectors will be in the rows of c
  pcps::xsyev('V', 'U', n, &c.at(0), n, &eval.at(0), &work.at(0), lwork, &rwork.at(0), info);
  if (info != 0)
    throw pcps::Exception("pcps::xsyev failed with error code %i in pcps::hf::diagonalize_fock") % info;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   constructs the integral arrays used in hartree fock calculations
///
/// \param[in]     ham      object containing the hamiltonian information
/// \param[in]     n        number of orbitals
/// \param[out]    oei_a    on exit, the one body integrals for alpha, size n*n
/// \param[out]    oei_b    on exit, the one body integrals for beta , size n*n
/// \param[out]    tei      on exit, the two body integrals, which are assumed real, in 1122 order, size n*n*n*n
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::build_ints(const pcps::Hamiltonian<S> & ham,
                                                  const int n,
                                                  std::vector<S> & oei_a,
                                                  std::vector<S> & oei_b,
                                                  std::vector<S> & tei)
{

  // quantum chemistry hamiltonian
  if ( ham.is_qc() ) {

    // one electron integrals
    for (int p = 0; p < n; p++)
    for (int q = 0; q < n; q++)
      oei_a[p*n+q] = oei_b[p*n+q] = ham.get_oei(p,q);

    // two electron integrals
    for (int p = 0; p < n; p++)
    for (int q = 0; q < n; q++)
    for (int r = 0; r < n; r++)
    for (int s = 0; s < n; s++)
      tei[ n*n*n*p + n*n*q + n*r + s ] = ham.get_tei(p,q,r,s);

  // hamiltonian built from individual operators
  } else {

    // zero the integrals
    pcps::xscal(n*n, pcps::zero<S>(), &oei_a.at(0), 1);
    pcps::xscal(n*n, pcps::zero<S>(), &oei_b.at(0), 1);
    pcps::xscal(n*n*n*n, pcps::zero<S>(), &tei.at(0), 1);

    // declare an iterator for looping over the hamiltonian's operators
    typename std::vector<typename pcps::OpBase<S>::const_ptr>::const_iterator op;

    // extract integrals from the hamiltonian's operators
    for (op = ham.ops().begin(); op != ham.ops().end(); op++) {

      // raising/lowering operators
      if (typeid(**op) == typeid(pcps::OpRL<S>)) {

        // check that we do not have any ab contributions
        if ( (**op).sites().at(0) % 2 != (**op).sites().at(1) % 2 )
          throw pcps::Exception("RL operators with one alpha and one beta site cannot be used in pcps::hf::build_ints");

        // add the operator's contribution to the one electron integrals
        const int site1 = (**op).sites().at(0)/2;
        const int site2 = (**op).sites().at(1)/2;
        if ((**op).sites().at(0) % 2 == 0 )
          oei_a[ n * site1 + site2 ] += (**op).coeff();
        else
          oei_b[ n * site1 + site2 ] += (**op).coeff();

      // double number operators
      } else if (typeid(**op) == typeid(pcps::OpNN<S>)) {

        // we currently only support on-site two electron terms like in the Hubbard model
        const int s0 = (**op).sites().at(0);
        const int s1 = (**op).sites().at(1);
        if ( ! ( s0 % 2 == 0 && s1 == s0 + 1 ) &&
             ! ( s0 % 2 == 1 && s1 == s0 - 1 )    )
          throw pcps::Exception("NN operators must act on one spatial site in pcps::hf::build_ints");

        // add the operator's contribution to the two electron integrals
        const int s = (**op).sites().at(0) / 2;
        tei[ n*n*n*s + n*n*s + n*s + s ] += (**op).coeff();

      // other operators are not currently supported
      } else {
        throw pcps::Exception("operator type \"%s\" not supported by pcps::hf::build_ints") % pcps::type_lookup(typeid(**op));
      }

    }

  }

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the alpha contribution to the hartree fock energy
///
/// \param[in]     n        number of orbitals
/// \param[in]     oei      the one body integrals, size n*n
/// \param[in]     f        the matrix for the fock operator, size n*n
/// \param[in]     rdm      the one body reduced density matrix, size n*n
///
/// \return the alpha contribution to the hartree fock energy
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> S pcps::hf::alpha_energy(const int n,
                                                      const std::vector<S> & oei,
                                                      const std::vector<S> & f,
                                                      const std::vector<S> & rdm)
{

  // initialize the energy
  S retval = pcps::zero<S>();

  // compute the energy
  for (int p = 0; p < n; p++)
  for (int q = 0; q < n; q++)
    retval += 0.5 * ( oei[p*n+q] + f[p*n+q] ) * rdm[p*n+q];

  // return the energy
  return retval;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   print the hartree fock orbitals and orbitals energies
///
/// \param[in]     n        number of orbitals
/// \param[in]     na       number of alpha particles
/// \param[in]     evals    the orbital energy eigenvalues, size n
/// \param[in]     c        the matrix of orbital coefficients,  size n*n
/// \param[in]     str      string for the particle type
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::print_orbitals(const int n,
                                                      const int na,
                                                      const std::vector<double> & evals,
                                                      const std::vector<S> & c,
                                                      const std::string & str)
{

  if (str.size() > 0)
    std::cout << boost::format("MO coefficients for %s:") % str << std::endl;
  else
    std::cout << boost::format("MO coefficients:") << std::endl;
  std::cout << std::endl;
  for (int i = 0; i < n; i++) {
    std::cout << boost::format("%20.9f  ") % evals[i];
    if (i < na)
      std::cout << boost::format("%4s  ") % "occ";
    else
      std::cout << boost::format("%4s  ") % "vir";
    for (int s = 0; s < n; s++)
      std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", c[i*n+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

}

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   computes the unrestricted hartree fock energy and orbital coefficients
///
/// \param[in]     n        number of orbitals
/// \param[in]     na       number of alpha particles
/// \param[in]     nb       number of beta  particles
/// \param[in]     ham      object containing the hamiltonian information
/// \param[in]     bias     for each orbital, a bias towards alpha or beta, size n
/// \param[in]     thresh   energy convergence threshold
/// \param[in]     maxit    maximum number of self consistent iterations
/// \param[in,out] ca       On entry, the initial guess for the alpha orbital coefficients,  size n*n.
///                         On exit, the optimized matrix of alpha orbital coefficients.
/// \param[in,out] cb       On entry, the initial guess for the beta orbital coefficients,  size n*n.
///                         On exit, the optimized matrix of beta orbital coefficients.
/// \param[out]    energy   On exit, the unrestricted hartree fock energy.
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class S> void pcps::hf::uhf(const int n,
                                           const int na,
                                           const int nb,
                                           const pcps::Hamiltonian<S> & ham,
                                           const std::vector<double> & bias,
                                           const double bias_decay,
                                           const double thresh,
                                           const int maxit,
                                           std::vector<S> & ca,
                                           std::vector<S> & cb,
                                           S & energy)
{

  // print a greeting
  std::cout << boost::format("starting unrestricted hartree fock iterations") << std::endl;
  std::cout << std::endl;

  // allocate and initialize the integrals
  std::vector<S> oei_a(n*n, pcps::zero<S>());
  std::vector<S> oei_b(n*n, pcps::zero<S>());
  std::vector<S> tei(n*n*n*n, pcps::zero<S>());
  pcps::hf::build_ints(ham, n, oei_a, oei_b, tei);

  // allocate fock matrices
  std::vector<S> fa(n*n, pcps::zero<S>());
  std::vector<S> fb(n*n, pcps::zero<S>());

  // allocate one body reduced density matrices
  std::vector<S> ra(n*n, pcps::zero<S>());
  std::vector<S> rb(n*n, pcps::zero<S>());

  // allocate orbital energy vectors
  std::vector<double> eval_a(n, 0.0);
  std::vector<double> eval_b(n, 0.0);

  // initialize the decay coefficient for the bias
  double decay = 1.0;

  // iterate
  bool converged = false;
  for (int iter = 0; iter < maxit; iter++) {

    // build one body density matrices
    pcps::hf::build_rdm1(n, na, ca, ra);
    pcps::hf::build_rdm1(n, nb, cb, rb);

    // build fock matrices
    pcps::hf::build_fock(n, ra, rb, oei_a, tei, fa);
    pcps::hf::build_fock(n, rb, ra, oei_b, tei, fb);

    // compute energy
    const S last_energy = energy;
    energy = pcps::unity<S>() * ham.enuc();
    energy += pcps::hf::alpha_energy(n, oei_a, fa, ra);
    energy += pcps::hf::alpha_energy(n, oei_b, fb, rb);

    // print energy
    std::cout << boost::format("%4i   %s") % iter % pcps::format_number<S>("%20.12f", energy) << std::endl;

    // check for convergence
    converged = ( iter > 0 && std::abs( energy - last_energy ) < thresh );
    if (converged) break;

    // add the bias to the fock matrices
    for (int p = 0; p < n; p++) {
      fa[n*p+p] -= pcps::unity<S>() * decay * bias.at(p);
      fb[n*p+p] += pcps::unity<S>() * decay * bias.at(p);
    }

    // diagonalize fock matrices to find the new orbital coefficients
    pcps::hf::diagonalize_fock(n, na, fa, eval_a, ca);
    pcps::hf::diagonalize_fock(n, nb, fb, eval_b, cb);

    // reduce the bias based on the user specified decay rate
    decay *= bias_decay;

  }

  // print results
  std::cout << std::endl;
  if (converged)
    std::cout << boost::format("unrestricted hartree fock has converged") << std::endl;
  else
    std::cout << boost::format("!!! unrestricted hartree fock failed to converge !!!") << std::endl;
  std::cout << std::endl;

  // print the mo coefficients and orbital energies
  pcps::hf::print_orbitals(n, na, eval_a, ca, "alpha");
  pcps::hf::print_orbitals(n, nb, eval_b, cb, "beta");

}

// explicitly instantiate the function
template void pcps::hf::uhf(const int n,
                            const int na,
                            const int nb,
                            const pcps::Hamiltonian<double> & ham,
                            const std::vector<double> & bias,
                            const double bias_decay,
                            const double thresh,
                            const int maxit,
                            std::vector<double> & ca,
                            std::vector<double> & cb,
                            double & energy);
template void pcps::hf::uhf(const int n,
                            const int na,
                            const int nb,
                            const pcps::Hamiltonian<std::complex<double> > & ham,
                            const std::vector<double> & bias,
                            const double bias_decay,
                            const double thresh,
                            const int maxit,
                            std::vector<std::complex<double> > & ca,
                            std::vector<std::complex<double> > & cb,
                            std::complex<double> & energy);

//-------------------------------------------------------------------------------
// pcps::hartree_fock_spinless -- computes the 1-body reduced density matrix of
//                                the Hartree Fock wavefunction in the local
//                                orbital basis.
//                                Note that the atomic orbitals are assumed to be
//                                orthonormal so that the overlap matrix is the
//                                identity matrix.
//-------------------------------------------------------------------------------

template <class S>
void pcps::hartree_fock_spinless(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham, std::vector<S> & rdm1) {

  // get the number of sites and particles
  const int nsites = userinp.nsites();
  const int nparticles = userinp.nparticles();

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // Start a timer
  pcps::Stopwatch sw; sw.start(); 

  // non-root processes wait to receive rdm from root process, then exit
  if (myrank != 0) {
    rdm1 = std::vector<S>(nsites*nsites, pcps::zero<S>());
    comm.Bcast((void*)&rdm1[0], rdm1.size() * sizeof(S), MPI::CHAR, 0);
    return;
  }

  // Print some information
  std::cout << boost::format("computing the 1-body RDM for the spinless Hartree-Fock wavefunction") << std::endl;
  std::cout << std::endl;
  std::cout << boost::format("%4s   %14s") % "iter" % "energy" << std::endl;

  // Initialize the orbital coefficients.
  // The coefficients are indexed by (mo_index * nsites + ao_index).
  std::vector<S> c(nsites*nsites, pcps::zero<S>());

  // Allocate the Fock matrix
  std::vector<S> f(nsites*nsites, pcps::zero<S>());

  // Initialize the 1-body RDM
  rdm1 = std::vector<S>(nsites*nsites, pcps::zero<S>());

  //// if restarting, load old rdm
  //if (userinp.restart()) {
  //  std::string file_path = userinp.workdir() + "rdm.bin";
  //  std::ifstream infile(file_path.c_str(), std::ios::in | std::ios::binary);
  //  if (!infile.is_open())
  //    throw pcps::Exception( (boost::format("failed to open file \"%s\" for reading old 1-body RDM") % file_path).str() );
  //  infile.read( (char *)(&rdm1[0]), rdm1.size() * sizeof(double) );
  //  infile.close();
  //}

  // Allocate arrays for use in diagonalizing the Fock matrix
  int info;
  const int lwork = nsites*10;
  std::vector<S> work(lwork);
  std::vector<double> rwork(3*nsites);
  std::vector<double> f_evals(nsites);

  //// Reorganize the two-electron integrals for use in matrix multiplication
  //std::vector<double> two_int(nsites*nsites*nsites*nsites, 0.0);
  //for (int s = 0; s < nsites; s++)
  //for (int r = 0; r < nsites; r++)
  //for (int q = 0; q < nsites; q++)
  //for (int p = 0; p < nsites; p++)
  //  two_int[nsites*(nsites*(nsites*s+r)+q)+p] = ham.get_te_coeff(p,r,s,q) - ham.get_te_coeff(p,q,s,r);

  // convert the hamiltonian's operators into 1- and 2-body integral tensors
  std::vector<S> one_int(nsites*nsites, pcps::zero<S>());
  std::vector<S> two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
  {
    std::vector<S> temp_one_int(nsites*nsites, pcps::zero<S>());
    std::vector<S> temp_two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
    for (typename std::vector<typename pcps::OpBase<S>::const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++)
      (**op).get_uhf_ints(nsites, temp_one_int, temp_two_int);
    for (int i = 0; i < temp_one_int.size(); i++)
      one_int[i] = temp_one_int[i];
    for (int i = 0; i < temp_two_int.size(); i++)
      two_int[i] = temp_two_int[i];
  }

  // Iteratively update the coefficients
  S energy = 1.0e100 * pcps::unity<S>();
  bool converged = false;
  for (int iter = 0; !converged; iter++) {

    // Build the Fock matrix
    const int n2 = nsites*nsites;
    pcps::xgemm('N', 'N', n2, 1, n2, 2.0 * pcps::unity<S>(), &two_int[0], n2, &rdm1[0], n2, pcps::zero<S>(), &f[0], n2);
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
      f[p*nsites+q] += one_int[p*nsites+q];

    // Compute the energy
    const S last_energy = energy;
    energy = pcps::zero<S>();
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
      energy = energy + 0.5 * rdm1[p*nsites+q] * ( f[p*nsites+q] + one_int[p*nsites+q] );

    // Print the energy for this iteration
    std::cout << boost::format("%4i   %20.12f") % iter % energy << std::endl;

    // Check for convergence
    converged = (std::abs(last_energy - energy) < 1.0e-9);
    if (converged || iter >= userinp.max_hf_iter()) break;

    // Diagonalize the Fock matrix to find the new orbital coefficients
    for (int i = 0; i < f.size(); i++)
      c[i] = f[i];
    pcps::xsyev('V', 'U', nsites, &c[0], nsites, &f_evals[0], &work[0], lwork, &rwork[0], info);
    if (info != 0)
      throw pcps::Exception( (boost::format("xsyev for spinless hartree fock failed with error code %i") % info).str() );

    // Compute the new 1-body RDM
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++) {
      rdm1[p*nsites+q] = pcps::zero<S>();
      for (int i = 0; i < nparticles; i++)
        rdm1[p*nsites+q] += pcps::conj( c[i*nsites+p] ) * c[i*nsites+q];
    }

  }

  // Print the result of the iteration
  if (converged)
    std::cout << std::endl
              << boost::format("spinless Hartree-Fock iterations converged!") << std::endl
              << std::endl;
  else
    std::cout << std::endl
              << boost::format("no convergence...") << std::endl
              << std::endl;

  // Print the solution's 1-body RDM in the AO basis
  std::cout << boost::format("1-body RDM:") << std::endl;
  std::cout << std::endl;
  for (int r = 0; r < nsites; r++) {
    for (int s = 0; s < nsites; s++)
      std::cout << pcps::format_number("%12.5f", rdm1[r*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

  // Print the 1-body RDM's trace
  S trace = pcps::zero<S>();
  for (int r = 0; r < nsites; r++)
    trace += rdm1[r*nsites+r];
  std::cout << "trace = " << pcps::format_number("%.5f", trace) << std::endl;
  std::cout << std::endl;

  // Print the solution's MO coefficients
  std::cout << boost::format("MO coefficients:") << std::endl;
  std::cout << std::endl;
  for (int i = 0; i < nsites; i++) {
    std::cout << boost::format("%20.9f  ") % f_evals[i];
    if (i < nparticles)
      std::cout << boost::format("%4s  ") % "occ";
    else
      std::cout << boost::format("%4s  ") % "vir";
    for (int s = 0; s < nsites; s++)
      std::cout << pcps::format_number("%12.5f", c[i*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

  // save 1-body RDM to disk
  {
    std::string file_path = userinp.workdir() + "rdm.bin";
    std::ofstream outfile(file_path.c_str(), std::ios::out | std::ios::binary);
    if (!outfile.is_open())
      throw pcps::Exception( (boost::format("failed to open file \"%s\" for saving 1-body RDM") % file_path).str() );
    outfile.write( (const char *)(&rdm1[0]), rdm1.size() * sizeof(S) );
    outfile.close();
  }

  // send the rdm to other processes
  comm.Bcast((void*)&rdm1[0], rdm1.size() * sizeof(S), MPI::CHAR, 0);
  
  // Print the elapsed time
  std::cout << boost::format("time to compute 1-body RDM for spinless Hartree-Fock: %.4f sec") % sw.elapsed_seconds() << std::endl;
  std::cout << std::endl;

}

// explicitly instantiate the function
template void pcps::hartree_fock_spinless(const pcps::Input &, const pcps::Hamiltonian< double               > &, std::vector< double               > &);
template void pcps::hartree_fock_spinless(const pcps::Input &, const pcps::Hamiltonian< std::complex<double> > &, std::vector< std::complex<double> > &);

//-------------------------------------------------------------------------------
// pcps::hartree_fock_restricted -- computes the 1-body reduced density matrix of
//                                  the Hartree Fock wavefunction in the local
//                                  orbital basis.
//                                  Note that the atomic orbitals are assumed to
//                                  be orthonormal so that the overlap matrix is
//                                  the identity matrix.
//                                  The output RDM is the full density matrix,
//                                  including both alpha and beta orbitals.
//                                  It therefore contains many zeros.
//-------------------------------------------------------------------------------

template <class S>
void pcps::hartree_fock_restricted(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham, std::vector<S> & c) {

  // check that the total number of sites and particles is even
  if (userinp.nsites() % 2 != 0)
    throw pcps::Exception( (boost::format("hartree_fock_restricted assumes there are an even number of lattice sites")).str() );
  if (userinp.nparticles() % 2 != 0)
    throw pcps::Exception( (boost::format("hartree_fock_restricted assumes there are an even number of particles")).str() );

  // find the number of alpha sites and particles
  const int nsites = userinp.nsites() / 2;
  const int nparticles = userinp.nparticles() / 2;

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // ensure the orbital coefficients are the correct size
  if (c.size() != nsites*nsites)
    throw pcps::Exception( (boost::format("wrong size for restricted orbital coefficient array in pcps::hartree_fock_restricted")).str() );

  // Start a timer
  pcps::Stopwatch sw; sw.start(); 

  // Print some information
  std::cout << boost::format("computing restricted orbital coefficients via a restricted Hartree-Fock calculation") << std::endl;
  std::cout << std::endl;
  std::cout << boost::format("%4s   %14s") % "iter" % "energy" << std::endl;

  // The coefficients are indexed by (mo_index * nsites + ao_index).

  // Allocate the Fock matrix
  std::vector<S> f(nsites*nsites, pcps::zero<S>());

  // Initialize the 1-body RDM
  std::vector<S> rdm1_alpha(nsites*nsites, pcps::zero<S>());
  for (int p = 0; p < nsites; p++)
  for (int q = 0; q < nsites; q++) {
    rdm1_alpha[p*nsites+q] = pcps::zero<S>();
    for (int i = 0; i < nparticles; i++)
      rdm1_alpha[p*nsites+q] += pcps::conj( c[i*nsites+p] ) * c[i*nsites+q];
  }

  // Allocate arrays for use in diagonalizing the Fock matrix
  int info;
  const int lwork = nsites*10;
  std::vector<S> work(lwork);
  std::vector<double> rwork(3*nsites);
  std::vector<double> f_evals(nsites);

  // convert the hamiltonian's operators into 1- and 2-body integral tensors
  std::vector<S> one_int(nsites*nsites, pcps::zero<S>());
  std::vector<S> two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
  if ( ham.is_qc() ) {
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
      one_int[nsites*i+j] = ham.get_oei(i,j);
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
    for (int k = 0; k < nsites; k++)
    for (int l = 0; l < nsites; l++)
      two_int[nsites*nsites*nsites*i + nsites*nsites*j + nsites*k + l] = ham.get_tei(i,j,k,l);
  } else {
    std::vector<S> temp_one_int(nsites*nsites, pcps::zero<S>());
    std::vector<S> temp_two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
    for (typename std::vector<typename pcps::OpBase<S>::const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++)
      (**op).get_rhf_ints(nsites, temp_one_int, temp_two_int);
    for (int i = 0; i < temp_one_int.size(); i++)
      one_int[i] = temp_one_int[i];
    for (int i = 0; i < temp_two_int.size(); i++)
      two_int[i] = temp_two_int[i];
  }

  // Iteratively update the coefficients
  S energy = 1.0e100 * pcps::unity<S>();
  bool converged = false;
  for (int iter = 0; !converged; iter++) {

    // Build the Fock matrix
    const int n2 = nsites*nsites;
    // coulomb term
    pcps::xgemm('N', 'N', n2, 1, n2, 2.0 * pcps::unity<S>(), &two_int[0], n2, &rdm1_alpha[0], n2, pcps::zero<S>(), &f[0], n2);
    // exchange term
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
    for (int r = 0; r < nsites; r++)
    for (int s = 0; s < nsites; s++)
      f[p*nsites+q] -= rdm1_alpha[r*nsites+s] * two_int[nsites*nsites*nsites*p + nsites*nsites*r + nsites*q + s];
    // one electron term
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
      f[p*nsites+q] += one_int[p*nsites+q];

    // Compute the energy
    const S last_energy = energy;
    energy = pcps::unity<S>() * ham.enuc();
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
      energy = energy + rdm1_alpha[p*nsites+q] * ( f[p*nsites+q] + one_int[p*nsites+q] );

    // Print the energy for this iteration
    std::cout << boost::format("%4i   %s") % iter % pcps::format_number<S>("%20.12f", energy) << std::endl;

    // Check for convergence
    converged = (std::abs(last_energy - energy) < 1.0e-9);
    if (converged || iter >= userinp.max_hf_iter()) break;

    // Diagonalize the Fock matrix to find the new orbital coefficients
    for (int i = 0; i < f.size(); i++)
      c[i] = f[i];
    pcps::xsyev('V', 'U', nsites, &c[0], nsites, &f_evals[0], &work[0], lwork, &rwork[0], info);
    if (info != 0)
      throw pcps::Exception( (boost::format("pcps::xsyev for restricted hartree fock failed with error code %i") % info).str() );

    //// tinker with the orbital coefficients when HOMO/LUMO have the same energy
    //if ( typeid(S) == typeid(std::complex<double>) && std::fabs(f_evals[nparticles-1] - f_evals[nparticles]) < 1.0e-9 )
    //  for (int p = 0; p < nsites; p++) {
    //    const S t1 = c[(nparticles-1)*nsites+p];
    //    const S t2 = c[(nparticles  )*nsites+p];
    //    const S i = std::sqrt(-1.0 * pcps::unity<S>());
    //    c[(nparticles-1)*nsites+p] = ( t1 + i * t2 ) / std::sqrt( 2.0 * pcps::unity<S>() );
    //    c[(nparticles  )*nsites+p] = ( t1 - i * t2 ) / std::sqrt( 2.0 * pcps::unity<S>() );
    //  }

    //// Print the new orbital coefficients
    //std::cout << std::endl;
    //std::cout << boost::format("MO coefficients:") << std::endl;
    //std::cout << std::endl;
    //for (int i = 0; i < nsites; i++) {
    //  std::cout << boost::format("%20.9f  ") % f_evals[i];
    //  if (i < nparticles)
    //    std::cout << boost::format("%4s  ") % "occ";
    //  else
    //    std::cout << boost::format("%4s  ") % "vir";
    //  for (int s = 0; s < nsites; s++)
    //    std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", c[i*nsites+s]);
    //  std::cout << std::endl;
    //}
    //std::cout << std::endl;

    // Compute the new 1-body RDM
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++) {
      rdm1_alpha[p*nsites+q] = pcps::zero<S>();
      for (int i = 0; i < nparticles; i++)
        rdm1_alpha[p*nsites+q] += pcps::conj( c[i*nsites+p] ) * c[i*nsites+q];
    }

    //// Print 1-body RDM in the AO basis
    //std::cout << boost::format("1-body alpha RDM:") << std::endl;
    //std::cout << std::endl;
    //for (int r = 0; r < nsites; r++) {
    //  for (int s = 0; s < nsites; s++)
    //    std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", rdm1_alpha[r*nsites+s]);
    //  std::cout << std::endl;
    //}
    //std::cout << std::endl;

  }

  // Print the result of the iteration
  if (converged)
    std::cout << std::endl
              << boost::format("restricted Hartree-Fock iterations converged!") << std::endl
              << std::endl;
  else
    std::cout << std::endl
              << boost::format("no convergence...") << std::endl
              << std::endl;

  // Print the solution's 1-body RDM in the AO basis
  std::cout << boost::format("1-body alpha RDM:") << std::endl;
  std::cout << std::endl;
  for (int r = 0; r < nsites; r++) {
    for (int s = 0; s < nsites; s++)
      std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", rdm1_alpha[r*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

  // Print the 1-body RDM's trace
  S trace = pcps::zero<S>();
  for (int r = 0; r < nsites; r++)
    trace += 2.0 * rdm1_alpha[r*nsites+r];
  std::cout << boost::format("trace = %s") % pcps::format_number<S>("%.5f", trace) << std::endl;
  std::cout << std::endl;

  // Print the solution's MO coefficients
  std::cout << boost::format("MO coefficients:") << std::endl;
  std::cout << std::endl;
  for (int i = 0; i < nsites; i++) {
    std::cout << boost::format("%20.9f  ") % f_evals[i];
    if (i < nparticles)
      std::cout << boost::format("%4s  ") % "occ";
    else
      std::cout << boost::format("%4s  ") % "vir";
    for (int s = 0; s < nsites; s++)
      std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", c[i*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

//  // convert the 1-body RDM back to the alpha and beta basis
//  rdm1 = std::vector<S>(4*nsites*nsites, pcps::zero<S>());
//  for (int r = 0; r < nsites; r++)
//  for (int s = 0; s < nsites; s++) {
//    rdm1[(2*r+0) * (2*nsites) + (2*s+0)] = rdm1_alpha[r*nsites+s];
//    rdm1[(2*r+1) * (2*nsites) + (2*s+1)] = rdm1_alpha[r*nsites+s];
//  }

  // Print the elapsed time
  std::cout << boost::format("restricted Hartree-Fock calculation took %.4f sec") % sw.elapsed_seconds() << std::endl;
  std::cout << std::endl;

}

// explicitly instantiate the function
template void pcps::hartree_fock_restricted(const pcps::Input &, const pcps::Hamiltonian< double               > &, std::vector< double               > &);
template void pcps::hartree_fock_restricted(const pcps::Input &, const pcps::Hamiltonian< std::complex<double> > &, std::vector< std::complex<double> > &);

//-------------------------------------------------------------------------------
// pcps::hartree_fock_restricted_quad -- computes the 1-body reduced density matrix of
//                                       the Hartree Fock wavefunction in the local
//                                       orbital basis.
//                                       Note that the atomic orbitals are assumed to
//                                       be orthonormal so that the overlap matrix is
//                                       the identity matrix.
//                                       The output RDM is the full density matrix,
//                                       including both alpha and beta orbitals.
//                                       It therefore contains many zeros.
//                                       This function uses a quadradically
//                                       convergent Newton-Raphson minimization
//                                       of the RHF energy.
//-------------------------------------------------------------------------------

template <class S>
void pcps::hartree_fock_restricted_quad(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham, std::vector<S> & rdm1) {

  if (typeid(S) != typeid(double))
    throw pcps::Exception( (boost::format("hartree_fock_restricted_quad only implemented for S = double")).str() );

  // check that the total number of sites and particles is even
  if (userinp.nsites() % 2 != 0)
    throw pcps::Exception( (boost::format("hartree_fock_restricted_quad assumes there are an even number of lattice sites")).str() );
  if (userinp.nparticles() % 2 != 0)
    throw pcps::Exception( (boost::format("hartree_fock_restricted_quad assumes there are an even number of particles")).str() );

  // find the number of alpha sites and particles
  const int nsites = userinp.nsites() / 2;
  const int nparticles = userinp.nparticles() / 2;

  // create a function object to orthonormalize a set of vectors via a graham schmidt procedure
  struct Local_Graham_Schmidt {
    // the vectors are assumed to be stored in the rows of the vecs matrix
    bool operator()(const int dim, const int nvec, S * const vecs) {
      bool independent = true;
      for (int i = 0; i < nvec; i++) {
        for (int j = 0; j < i; j++) {
          S inner_prod = pcps::zero<S>();
          for (int p = 0; p < dim; p++)
            inner_prod += pcps::conj(vecs[j*dim+p]) * vecs[i*dim+p];
          pcps::xaxpy(dim, -inner_prod, &vecs[j*dim], 1, &vecs[i*dim], 1);
        }
        S norm = pcps::zero<S>();
        for (int p = 0; p < dim; p++)
          norm += pcps::conj(vecs[i*dim+p]) * vecs[i*dim+p];
        norm = std::sqrt(norm);
        if (std::abs(norm) < 1.0e-10) {
          independent = false;
          break;
        }
        pcps::xscal(dim, pcps::unity<S>() / norm, &vecs[i*dim], 1);
      }
      return independent;
    }
  };
  Local_Graham_Schmidt local_graham_schmidt;

  // create a function object that maps a pair of integers i,j (0 <= i <= j < n) to a compound index in the range 0, 1, ..., n*(n+1)/2-1
  class Local_Compound {
    private:
      const int _n;
    public:
      Local_Compound(const int n) : _n(n) {}
      int operator()(const int i, const int j) const {
        assert(i <= j);
        return ( 2 * _n - i + 1 ) * i / 2 + j - i;
      }
  };
  const Local_Compound local_cmpd(nparticles);

  // create a function object to compute the alpha rdm given a set of coefficients
  struct Local_Build_RDM {
    void operator()(const int nsites, const int nparticles, const S * coeff, S * rdm) {
      pcps::xscal(nsites*nsites, pcps::zero<S>(), &rdm[0], 1);
      for (int p = 0; p < nsites; p++)
      for (int q = 0; q < nsites; q++)
      for (int i = 0; i < nparticles; i++)
        rdm[p*nsites+q] += pcps::conj( coeff[i*nsites+p] ) * coeff[i*nsites+q];
    }
  };
  Local_Build_RDM local_build_rdm;

  // create a function object to compute the RHF fock matrix given a set of coefficients and integrals
  struct Local_Build_RHF_Fock_Matrix {
    void operator()(const int nsites, const int nparticles, const S * rdm, const S * oei, const S * tei, S * fmat) {
      const int n2 = nsites*nsites;
      pcps::xgemm('N', 'N', n2, 1, n2, 2.0 * pcps::unity<S>(), &tei[0], n2, &rdm[0], n2, pcps::zero<S>(), &fmat[0], n2);
      for (int p = 0; p < nsites; p++)
      for (int q = 0; q < nsites; q++)
      for (int r = 0; r < nsites; r++)
      for (int s = 0; s < nsites; s++)
        fmat[p*nsites+q] -= rdm[r*nsites+s] * tei[nsites*(nsites*(nsites*s+q)+r)+p];
      for (int p = 0; p < nsites; p++)
      for (int q = 0; q < nsites; q++)
        fmat[p*nsites+q] += oei[p*nsites+q];
    }
  };
  Local_Build_RHF_Fock_Matrix local_build_rhf_fock_matrix;

  // create a function object to compute the RHF energy given a set of coefficients and integrals
  struct Local_RHF_Energy {
    // requires work array of size 2*nsites*nsites
    S operator()(const int nsites, const int nparticles, const S * rdm, const S * fock_matrix, const S * oei) {
      // compute energy
      S retval = pcps::zero<S>();
      for (int p = 0; p < nsites; p++)
      for (int q = 0; q < nsites; q++)
        retval += rdm[p*nsites+q] * ( oei[p*nsites+q] + fock_matrix[p*nsites+q] );
      return retval;
    }
  };
  Local_RHF_Energy local_rhf_energy;

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // non-root processes wait to receive rdm from root process, then exit
  if (myrank != 0) {
    rdm1 = std::vector<S>(4*nsites*nsites, pcps::zero<S>());
    comm.Bcast((void*)&rdm1[0], rdm1.size() * sizeof(S), MPI::CHAR, 0);
    return;
  }

  // Start a timer
  pcps::Stopwatch sw; sw.start(); 

  // Print some information
  std::cout << boost::format("computing the 1-body RDM for the restricted Hartree-Fock wavefunction via Newton-Raphson") << std::endl;
  std::cout << std::endl;

  // Initialize the orbital coefficients.
  // The coefficients are indexed by (mo_index * nsites + ao_index).
  std::vector<S> c(nsites*nsites, pcps::zero<S>());

  // Initialize the 1-body RDM
  std::vector<S> rdm1_alpha(nsites*nsites, pcps::zero<S>());

  // Compute the Hessian's dimension
  const int hess_dim =  nsites * nparticles + nparticles * (nparticles+1) / 2;

  // Allocate the Fock matrix
  std::vector<S> f(nsites*nsites, pcps::zero<S>());

  // Allocate the Hesian matrix
  std::vector<S> hess(hess_dim*hess_dim, pcps::zero<S>());

  // Allocate the residual and update vectors
  std::vector<S> res(hess_dim, pcps::zero<S>());
  std::vector<S> update(hess_dim, pcps::zero<S>());

  // Allocate a vector for the lagrange multipliers
  std::vector<S> mult(nparticles * (nparticles+1) / 2, pcps::zero<S>());

  // Allocate arrays for use in matrix diagonalization
  int info;
  const int lwork = std::max( std::max(10*hess_dim, 10*nsites*nsites), nsites*nsites*nsites*nsites );
  std::vector<S> work(lwork, pcps::zero<S>());
  std::vector<double> rwork(3*hess_dim, 0.0);
  std::vector<double> evals(hess_dim, 0.0);

  // convert the hamiltonian's operators into 1- and 2-body integral tensors
  std::vector<S> one_int(nsites*nsites, pcps::zero<S>());
  std::vector<S> two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
  {
    std::vector<S> temp_one_int(nsites*nsites, pcps::zero<S>());
    std::vector<S> temp_two_int(nsites*nsites*nsites*nsites, pcps::zero<S>());
    for (typename std::vector<typename pcps::OpBase<S>::const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++)
      (**op).get_rhf_ints(nsites, temp_one_int, temp_two_int);
    for (int i = 0; i < temp_one_int.size(); i++)
      one_int[i] = temp_one_int[i];
    for (int i = 0; i < temp_two_int.size(); i++)
      two_int[i] = temp_two_int[i];
  }

  // Allocate arrays for the MO integrals and integral transformation intermediates
  std::vector<S> one_int_mo(nsites*nsites, pcps::zero<S>());
  std::vector<S> two_int_mo(nsites*nsites*nsites*nsites, pcps::zero<S>());
  S * two_int_xx = &work[0];

  // use the free particle solution as the initial guess
  {

    // zero the density matrix
    pcps::xscal(nsites*nsites, pcps::zero<S>(), &rdm1_alpha[0], 1);

    // add random noise to the density matrix
    //for (int p = 0; p < nsites; p++)
    //for (int q = p; q < nsites; q++) {
    //  rdm1_alpha[p*nsites+q] += 0.5 * pcps::random_number<S>();
    //  rdm1_alpha[q*nsites+p] = pcps::conj(rdm1_alpha[p*nsites+q]);
    //}

    // Build the Fock matrix
    local_build_rhf_fock_matrix(nsites, nparticles, &rdm1_alpha[0], &one_int[0], &two_int[0], &c[0]);

    // Diagonalize the Fock matrix to find the orbital coefficients
    pcps::xsyev('V', 'U', nsites, &c[0], nsites, &evals[0], &work[0], lwork, &rwork[0], info);
    if (info != 0)
      throw pcps::Exception( (boost::format("pcps::xsyev for restricted hartree fock failed with error code %i") % info).str() );

    // use the fock matrix eigenvalues as an initial guess for the lagrange multipliers
    for (int i = 0; i < nparticles; i++)
      mult[local_cmpd(i,i)] = -2.0 * pcps::unity<S>() * evals[i];

  }

  bool converged = false;
  for (int super_iter = 0; super_iter < 4; super_iter++) {

    // iterate the newton raphson method
    converged = false;
    for (int iter = 0; iter < userinp.max_hf_iter(); iter++) {

      // Compute the 1-body RDM
      local_build_rdm(nsites, nparticles, &c[0], &rdm1_alpha[0]);

      // Build the Fock matrix
      local_build_rhf_fock_matrix(nsites, nparticles, &rdm1_alpha[0], &one_int[0], &two_int[0], &f[0]);

      // Compute the energy (actually the value of the Lagrangian, which is equal to the energy at its minimum)
      const S energy = local_rhf_energy(nsites, nparticles, &rdm1_alpha[0], &f[0], &one_int[0]);

      // Compute the Lagrangian's derivatives, which we refer to as the residuals
      pcps::xscal(hess_dim, pcps::zero<S>(), &res[0], 1);
      for (int i = 0; i < nparticles; i++)
      for (int p = 0; p < nsites; p++) {
        S & value = res[i*nsites+p];
        for (int q = 0; q < nsites; q++)
          value += 4.0 * c[i*nsites+q] * f[p*nsites+q];
        for (int j = i; j < nparticles; j++)
          value += mult[local_cmpd(i, j)] * c[j*nsites+p];
        for (int j = 0; j <= i; j++)
          value += mult[local_cmpd(j, i)] * c[j*nsites+p];
      }
      for (int j = 0; j < nparticles; j++)
      for (int i = 0; i <= j; i++) {
        S & value = res[nsites*nparticles + local_cmpd(i,j)];
        value -= pcps::unity<S>() * (i == j ? 1.0 : 0.0);
        for (int p = 0; p < nsites; p++)
          value += pcps::conj(c[i*nsites+p]) * c[j*nsites+p];
      }

      // Compute the Lagrangian's value
      S lagrangian = energy;
      for (int j = 0; j < nparticles; j++)
      for (int i = 0; i <= j; i++)
        lagrangian += mult[local_cmpd(i,j)] * res[nsites*nparticles + local_cmpd(i,j)];

      //// Print the coefficient and lagrange multiplier values
      //std::cout << "values:" << std::endl;
      //for (int i = 0; i < nsites*nparticles; i++)
      //  std::cout << boost::format("%20.12e") % c[i] << std::endl;
      //for (int i = 0; i < nparticles*(nparticles+1)/2; i++)
      //  std::cout << boost::format("%20.12e") % mult[i] << std::endl;
      //std::cout << std::endl;

      //// Print the residual
      //std::cout << "residual:" << std::endl;
      //for (int i = 0; i < hess_dim; i++)
      //  std::cout << boost::format("%20.12e") % res[i] << std::endl;
      //std::cout << std::endl;

      // check convergence
      S res_norm = pcps::zero<S>();
      for (int i = 0; i < hess_dim; i++)
        res_norm += pcps::square_norm(res[i]);
      res_norm = std::sqrt(res_norm);
      converged = ( std::abs(res_norm) < 1.0e-9 );

      // Print the iteration information
      std::cout << boost::format("iteration %4i        lagrangian = %s       residual norm = %.2e")
                     % iter
                     % pcps::format_number("%20.12f", lagrangian)
                     % std::abs(res_norm)
                << std::endl;
      //std::cout << std::endl;

      // end iteration if converged
      if (converged) break;

      // Compute the Hessian matrix, which contains the second derivatives of the Lagrangian
      pcps::xscal(hess_dim*hess_dim, pcps::zero<S>(), &hess[0], 1);
      for (int J = 0; J < nparticles; J++)
      for (int Q = 0; Q < nsites; Q++)
      for (int I = 0; I <= J; I++)
      for (int P = 0; P < nsites; P++) {
        const int IP = I*nsites+P;
        const int JQ = J*nsites+Q;
        const int cmpd = IP*hess_dim+JQ;
        S & value = hess[cmpd];
        if (P == Q)
          value += mult[local_cmpd(I,J)];
        if (I == J) {
          value += 4.0 * one_int[P*nsites+Q];
          for (int r = 0; r < nsites; r++)
          for (int s = 0; s < nsites; s++)
            value += rdm1_alpha[r*nsites+s] * (   8.0 * two_int[nsites*(nsites*(nsites*s+r)+Q)+P]
                                                - 4.0 * two_int[nsites*(nsites*(nsites*s+Q)+r)+P] );
        }
        for (int q = 0; q < nsites; q++)
        for (int s = 0; s < nsites; s++)
          value += c[I*nsites+q] * c[J*nsites+s] * (   16.0 * two_int[nsites*(nsites*(nsites*s+Q)+q)+P]
                                                     -  4.0 * two_int[nsites*(nsites*(nsites*s+q)+Q)+P]
                                                     -  4.0 * two_int[nsites*(nsites*(nsites*Q+q)+s)+P] );
        hess[JQ*hess_dim+IP] = value;
      }
      for (int J = 0; J < nparticles; J++)
      for (int I = 0; I <= J; I++) {
        for (int Q = 0; Q < nsites; Q++) {
          {
            const int K = I;
            const int KQ = K*nsites+Q;
            hess[KQ*hess_dim + (nsites*nparticles + local_cmpd(I,J))] += c[J*nsites+Q];
            hess[(nsites*nparticles + local_cmpd(I,J))*hess_dim + KQ] += c[J*nsites+Q];
          }
          {
            const int K = J;
            const int KQ = K*nsites+Q;
            hess[KQ*hess_dim + (nsites*nparticles + local_cmpd(I,J))] += c[I*nsites+Q];
            hess[(nsites*nparticles + local_cmpd(I,J))*hess_dim + KQ] += c[I*nsites+Q];
          }
        }
      }

      // Diagonalize the Hessian (eigenvectors will be stored in the rows of hess)
      pcps::xsyev('V', 'U', hess_dim, &hess[0], hess_dim, &evals[0], &work[0], lwork, &rwork[0], info);
      if (info != 0)
        throw pcps::Exception( (boost::format("pcps::xsyev for restricted hartree fock failed with error code %i") % info).str() );

      //// Print the Hessian's eigenvalues
      //std::cout << "hessian eigenvalues:" << std::endl;
      //for (int i = 0; i < hess_dim; i++)
      //  std::cout << boost::format("%20.12e") % evals[i] << std::endl;
      //std::cout << std::endl;

      // Compute the new coefficient and lagrange multiplier values
      pcps::xgemm('C', 'N', hess_dim, 1, hess_dim, - pcps::unity<S>(), &hess[0], hess_dim, &res[0], hess_dim,
                  pcps::zero<S>(), &update[0], hess_dim);
      for (int i = 0; i < hess_dim; i++)
        update[i] /= evals[i];
      pcps::xgemm('N', 'N', hess_dim, 1, hess_dim, pcps::unity<S>(), &hess[0], hess_dim, &update[0], hess_dim,
                  pcps::zero<S>(), &res[0], hess_dim);
      for (int i = 0; i < hess_dim; i++)
        update[i] = res[i];
      const double step_size = 1.0;
      for (int i = 0; i < nparticles; i++)
      for (int p = 0; p < nsites; p++)
        c[i*nsites+p] += step_size * update[i*nsites+p];
      for (int j = 0; j < nparticles; j++)
      for (int i = 0; i <= j; i++)
        mult[local_cmpd(i,j)] += step_size * update[nsites*nparticles + local_cmpd(i,j)];

    }
    std::cout << std::endl;

    // Compute the 1-body RDM
    local_build_rdm(nsites, nparticles, &c[0], &rdm1_alpha[0]);

    //// Print the 1-body RDM in the AO basis
    //std::cout << boost::format("1-body alpha RDM:") << std::endl;
    //std::cout << std::endl;
    //for (int r = 0; r < nsites; r++) {
    //  for (int s = 0; s < nsites; s++)
    //    std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", rdm1_alpha[r*nsites+s]);
    //  std::cout << std::endl;
    //}
    //std::cout << std::endl;

    //// Print the occupied orbital coefficients
    //std::cout << boost::format("MO coefficients:") << std::endl;
    //std::cout << std::endl;
    //for (int i = 0; i < nparticles; i++) {
    //  std::cout << boost::format("%20s  ") % "";
    //  if (i < nparticles)
    //    std::cout << boost::format("%4s  ") % "occ";
    //  else
    //    std::cout << boost::format("%4s  ") % "vir";
    //  for (int s = 0; s < nsites; s++)
    //    std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", c[i*nsites+s]);
    //  std::cout << std::endl;
    //}
    //std::cout << std::endl;

    // Compute the virtual orbital coefficients (these are arbitrary but must be orthonormal wrt occupied orbitals and each other)
    while (true) {

      // assign arbitrary values to the virtual orbital coefficients
      for (int i = nparticles; i < nsites; i++)
      for (int p = 0; p < nsites; p++)
        c[i*nsites+p] = pcps::unity<S>() * pcps::random_number<double>();

      // enforce orthonormality through graham schmidt
      const bool linearly_independant = local_graham_schmidt(nsites, nsites, &c[0]);

      // If the vectors were linearly independent, we are finished.
      // Otherwise, choose a new set of random values for the virtual orbital coefficients
      if (linearly_independant) break;

    }

    // Compute integrals in MO basis
    pcps::xscal(nsites*nsites, pcps::zero<S>(), &one_int_mo[0], 1);
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
      one_int_mo[i*nsites+j] += c[i*nsites+p] * c[j*nsites+q] * one_int[p*nsites+q];
    pcps::xscal(nsites*nsites*nsites*nsites, pcps::zero<S>(), &two_int_xx[0], 1);
    for (int i = 0; i < nsites; i++)
    for (int p = 0; p < nsites; p++)
    for (int q = 0; q < nsites; q++)
    for (int r = 0; r < nsites; r++)
    for (int s = 0; s < nsites; s++)
      two_int_xx[nsites*(nsites*(nsites*s+r)+q)+i] += c[i*nsites+p] *    two_int[nsites*(nsites*(nsites*s+r)+q)+p];
    pcps::xscal(nsites*nsites*nsites*nsites, pcps::zero<S>(), &two_int_mo[0], 1);
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
    for (int q = 0; q < nsites; q++)
    for (int r = 0; r < nsites; r++)
    for (int s = 0; s < nsites; s++)
      two_int_mo[nsites*(nsites*(nsites*s+r)+j)+i] += c[j*nsites+q] * two_int_xx[nsites*(nsites*(nsites*s+r)+q)+i];
    pcps::xscal(nsites*nsites*nsites*nsites, pcps::zero<S>(), &two_int_xx[0], 1);
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
    for (int k = 0; k < nsites; k++)
    for (int r = 0; r < nsites; r++)
    for (int s = 0; s < nsites; s++)
      two_int_xx[nsites*(nsites*(nsites*s+k)+j)+i] += c[k*nsites+r] * two_int_mo[nsites*(nsites*(nsites*s+r)+j)+i];
    pcps::xscal(nsites*nsites*nsites*nsites, pcps::zero<S>(), &two_int_mo[0], 1);
    for (int i = 0; i < nsites; i++)
    for (int j = 0; j < nsites; j++)
    for (int k = 0; k < nsites; k++)
    for (int l = 0; l < nsites; l++)
    for (int s = 0; s < nsites; s++)
      two_int_mo[nsites*(nsites*(nsites*l+k)+j)+i] += c[l*nsites+s] * two_int_xx[nsites*(nsites*(nsites*s+k)+j)+i];

    // Compute energy using MO integrals
    S energy_check = pcps::zero<S>();
    for (int i = 0; i < nparticles; i++)
      energy_check += 2.0 * one_int_mo[i*nsites+i];
    for (int i = 0; i < nparticles; i++)
    for (int j = 0; j < nparticles; j++)
      energy_check += 2.0 * two_int_mo[nsites*(nsites*(nsites*j+j)+i)+i] - two_int_mo[nsites*(nsites*(nsites*j+i)+j)+i];

    // Print the energy
    std::cout << boost::format("energy computed via MO integrals = %s") % pcps::format_number("%20.12f", energy_check) << std::endl;
    std::cout << std::endl;

    // determine the dimension of the stability matrix
    const int stab_dim = nparticles * (nsites - nparticles);

    // construct the stability matrix
    pcps::xscal(stab_dim*stab_dim, pcps::zero<S>(), &hess[0], 1);
    for (int i = 0; i < nparticles; i++)
    for (int j = 0; j < nparticles; j++)
    for (int a = nparticles; a < nsites; a++)
    for (int b = nparticles; b < nsites; b++) {
      const int ia = i * (nsites-nparticles) + a - nparticles;
      const int jb = j * (nsites-nparticles) + b - nparticles;
      S & value = hess[ia * stab_dim + jb];
      value += 10.0 * two_int_mo[nsites*(nsites*(nsites*b+j)+a)+i] - 4.0 * two_int_mo[nsites*(nsites*(nsites*a+j)+b)+i];
      if        (a != b && i != j) {
        value += 2.0 * two_int_mo[nsites*(nsites*(nsites*b+j)+a)+i] - 2.0 * two_int_mo[nsites*(nsites*(nsites*b+a)+j)+i];
      } else if (a != b && i == j) {
        value += 2.0 * one_int_mo[a*nsites+b] - 2.0 * two_int_mo[nsites*(nsites*(nsites*j+i)+b)+a] + 2.0 * two_int_mo[nsites*(nsites*(nsites*b+j)+a)+i];
        for (int k = 0; k < nparticles; k++)
          value += 4.0 * two_int_mo[nsites*(nsites*(nsites*k+k)+b)+a] - 2.0 * two_int_mo[nsites*(nsites*(nsites*b+k)+a)+k];
      } else if (a == b && i != j) {
        value -= 2.0 * one_int_mo[i*nsites+j];
        value -= 2.0 * two_int_mo[nsites*(nsites*(nsites*b+a)+j)+i] - 2.0 * two_int_mo[nsites*(nsites*(nsites*b+j)+a)+i];
        for (int k = 0; k < nparticles; k++)
          value -= 4.0 * two_int_mo[nsites*(nsites*(nsites*k+k)+j)+i] - 2.0 * two_int_mo[nsites*(nsites*(nsites*k+j)+k)+i];
      } else if (a == b && i == j) {
        value -= 2.0 * energy_check;
        value += 2.0 * one_int_mo[a*nsites+a];
        value -= 2.0 * one_int_mo[i*nsites+i];
        for (int k = 0; k < nparticles; k++)
          value += 4.0 * one_int_mo[k*nsites+k];
        for (int k = 0; k < nparticles; k++)
        for (int l = 0; l < nparticles; l++) {
          const double factor = ( k == i || l == i ? 1.0 : 2.0 );
          value += factor * ( two_int_mo[nsites*(nsites*(nsites*l+l)+k)+k] - two_int_mo[nsites*(nsites*(nsites*l+k)+l)+k] );
        }
        for (int k = 0; k < nparticles; k++) {
          const double factor = ( k == i ? 1.0 : 2.0 );
          value += 2.0 * factor * two_int_mo[nsites*(nsites*(nsites*a+a)+k)+k] - 2.0 * (factor-1.0) * two_int_mo[nsites*(nsites*(nsites*a+k)+a)+k];
        }
        for (int k = 0; k < nparticles; k++)
        for (int l = 0; l < nparticles; l++) {
          if (k == i) continue;
          value += 2.0 * two_int_mo[nsites*(nsites*(nsites*l+l)+k)+k];
        }
      }
    }

    // Diagonalize the stability matrix (eigenvectors will be stored in the rows of hess)
    pcps::xsyev('V', 'U', stab_dim, &hess[0], stab_dim, &evals[0], &work[0], lwork, &rwork[0], info);
    if (info != 0)
      throw pcps::Exception( (boost::format("pcps::xsyev for restricted hartree fock failed with error code %i") % info).str() );

    // Print the stability matrix's eigenvalues
    std::cout << "stability matrix eigenvalues:" << std::endl;
    std::cout << std::endl;
    for (int i = 0; i < stab_dim; i++)
      std::cout << boost::format("%20.12e") % evals[i] << std::endl;
    std::cout << std::endl;

    // if the solution is unstable, change the coefficients to lower the energy
    if (evals[0] < 0.0) {

      // inform the user that the solution is unstable
      std::cout << boost::format("RHF solution is not stable.  Moving along an unstable eigenvector.") << std::endl;
      std::cout << std::endl;

      // try various step lengths to find a lower energy
      S last_energy = energy_check;
      S energy_after_step = pcps::zero<S>();
      double step_size;
      for (step_size = 0.01; step_size < 10.0; step_size += 0.01) {

        // take a step along the direction of the eigenvector with lowest eigenvalue
        pcps::xcopy(nsites*nsites, &c[0], 1, &work[0], 1);
        for (int i = 0; i < nparticles; i++)
        for (int a = nparticles; a < nsites; a++) {
          const int ia = i * (nsites-nparticles) + a - nparticles;
          pcps::xaxpy(nsites, step_size * hess[ia], &work[a*nsites], 1, &work[i*nsites], 1);
        }

        // enforce orthonormality through graham schmidt
        const bool linearly_independant = local_graham_schmidt(nsites, nsites, &work[0]);
        if (!linearly_independant)
          throw pcps::Exception("encountered a linear dependency while moving along an unstable eigenvector");

        // Compute the 1-body RDM
        local_build_rdm(nsites, nparticles, &work[0], &rdm1_alpha[0]);

        // Build the Fock matrix
        local_build_rhf_fock_matrix(nsites, nparticles, &rdm1_alpha[0], &one_int[0], &two_int[0], &f[0]);

        // Compute the energy
        energy_after_step = local_rhf_energy(nsites, nparticles, &rdm1_alpha[0], &f[0], &one_int[0]);

        // if the energy went up, stop iterating
        if (pcps::real(energy_after_step) > pcps::real(last_energy))
          break;

        // remember the energy for this step length
        last_energy = energy_after_step;

      }

      // save the new orbital coefficients
      pcps::xcopy(nsites*nsites, &work[0], 1, &c[0], 1);

      // Print the energy
      std::cout << boost::format("step size = %5.2f     new energy = %s") % step_size % pcps::format_number("%20.12f", energy_after_step) << std::endl;
      std::cout << std::endl;

    // if the solution was stable, stop iterating
    } else {
      break;
    }

  }

  // Print the result of the iteration
  if (converged)
    std::cout << boost::format("restricted Hartree-Fock iterations converged!") << std::endl
              << std::endl;
  else
    std::cout << boost::format("no convergence...") << std::endl
              << std::endl;

  // Print the solution's 1-body RDM in the AO basis
  std::cout << boost::format("1-body alpha RDM:") << std::endl;
  std::cout << std::endl;
  for (int r = 0; r < nsites; r++) {
    for (int s = 0; s < nsites; s++)
      std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", rdm1_alpha[r*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

  // Print the 1-body RDM's trace
  S trace = pcps::zero<S>();
  for (int r = 0; r < nsites; r++)
    trace += 2.0 * rdm1_alpha[r*nsites+r];
  std::cout << boost::format("trace = %s") % pcps::format_number<S>("%.5f", trace) << std::endl;
  std::cout << std::endl;

  // Print the solution's MO coefficients
  std::cout << boost::format("MO coefficients:") << std::endl;
  std::cout << std::endl;
  for (int i = 0; i < nsites; i++) {
    if (i < nparticles)
      std::cout << boost::format("%4s  ") % "occ";
    else
      std::cout << boost::format("%4s  ") % "vir";
    for (int s = 0; s < nsites; s++)
      std::cout << boost::format(" %s") % pcps::format_number<S>("%12.5f", c[i*nsites+s]);
    std::cout << std::endl;
  }
  std::cout << std::endl;

  // check that the orbitals are orthanormal
  for (int i = 0; i < nsites; i++)
  for (int j = 0; j < nsites; j++) {
    S ip = ( i == j ? pcps::unity<S>() : pcps::zero<S>() );
    for (int p = 0; p < nsites; p++)
      ip -= pcps::conj(c[i*nsites+p]) * c[j*nsites+p];
    if (std::abs(ip) > 1.0e-10)
      throw pcps::Exception( (boost::format("orbitals %i and %i failed orthonormality by %.2e") % i % j % std::abs(ip)).str() );
  }

  // convert the 1-body RDM back to the alpha and beta basis
  rdm1 = std::vector<S>(4*nsites*nsites, pcps::zero<S>());
  for (int r = 0; r < nsites; r++)
  for (int s = 0; s < nsites; s++) {
    rdm1[(2*r+0) * (2*nsites) + (2*s+0)] = rdm1_alpha[r*nsites+s];
    rdm1[(2*r+1) * (2*nsites) + (2*s+1)] = rdm1_alpha[r*nsites+s];
  }

  // send the rdm to other processes
  comm.Bcast((void*)&rdm1[0], rdm1.size() * sizeof(S), MPI::CHAR, 0);
   
  // Print the elapsed time
  std::cout << boost::format("time to compute 1-body RDM for restricted Hartree-Fock: %.4f sec") % sw.elapsed_seconds() << std::endl;
  std::cout << std::endl;

}

// explicitly instantiate the function
template void pcps::hartree_fock_restricted_quad(const pcps::Input &,
                                                 const pcps::Hamiltonian<double> &,
                                                 std::vector<double> &);
template void pcps::hartree_fock_restricted_quad(const pcps::Input &,
                                                 const pcps::Hamiltonian< std::complex<double> > &,
                                                 std::vector< std::complex<double> > &);

////-------------------------------------------------------------------------------
//// pcps::free_particle_ref -- computes the 1-body reduced density matrix of
////                            the solution to the free particle Hamiltonian,
////                            defined by setting the two-body integrals to zero
////                            Note that the atomic orbitals are assumed to be
////                            orthonormal so that the overlap matrix is the
////                            identity matrix.
////-------------------------------------------------------------------------------
//
//void pcps::free_particle_ref(const pcps::Input & userinp,
//                             const pcps::Hamiltonian & ham,
//                             std::vector<double> & rdm1)
//{
//  const int nsites = userinp.nsites();
//  const int nparticles = userinp.nparticles();
//
//  // Start a timer
//  pcps::Stopwatch sw; sw.start(); 
//
//  // Print some information
//  std::cout << boost::format("computing the 1-body RDM for the free particle wavefunction") << std::endl;
//  std::cout << std::endl;
//
//  // Initialize the orbital coefficients.
//  // The coefficients are indexed by (mo_index * nsites + ao_index).
//  std::vector<double> c(nsites*nsites, 0.0);
//
//  // Allocate the Fock matrix
//  std::vector<double> f(nsites*nsites, 0.0);
//
//  // Initialize the 1-body RDM
//  rdm1 = std::vector<double>(nsites*nsites, 0.0);
//
//  // convert the hamiltonian's operators into 1- and 2-body integral tensors
//  std::vector<double> one_int(nsites*nsites, 0.0);
//  std::vector<double> two_int(nsites*nsites*nsites*nsites, 0.0);
//  for (std::vector<pcps::OpBase_const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++)
//    (**op).get_uhf_ints(nsites, one_int, two_int);
//
//  // Allocate arrays for use in diagonalizing the Fock matrix
//  int info;
//  const int lwork = nsites*10;
//  std::vector<double> work(lwork);
//  std::vector<double> f_evals(nsites);
//
//  // Build the Fock matrix
//  for (int p = 0; p < nsites; p++)
//  for (int q = 0; q < nsites; q++)
//    f[p*nsites+q] += one_int[p*nsites+q];
//
//  // Diagonalize the Fock matrix to find the new orbital coefficients
//  for (int i = 0; i < f.size(); i++)
//    c[i] = f[i];
//  pcps::dsyev('V', 'U', nsites, &c[0], nsites, &f_evals[0], &work[0], lwork, info);
//  if (info != 0)
//    throw pcps::Exception( (boost::format("dsyev for free particle fock failed with error code %i") % info).str() );
//
//  // Compute the new 1-body RDM
//  pcps::dgemm('N', 'T', nsites, nsites, nparticles, 1.0, &c[0], nsites, &c[0], nsites, 0.0, &rdm1[0], nsites);
//
//  //// Reorganize the two-electron integrals for use in matrix multiplication
//  //std::vector<double> two_int(nsites*nsites*nsites*nsites, 0.0);
//  //for (int s = 0; s < nsites; s++)
//  //for (int r = 0; r < nsites; r++)
//  //for (int q = 0; q < nsites; q++)
//  //for (int p = 0; p < nsites; p++)
//  //  two_int[nsites*(nsites*(nsites*s+r)+q)+p] = ham.get_te_coeff(p,r,s,q) - ham.get_te_coeff(p,q,s,r);
//
//  // Compute the energy
//  const int n2 = nsites*nsites;
//  pcps::dgemm('N', 'N', n2, 1, n2, 1.0, &two_int[0], n2, &rdm1[0], n2, 0.0, &f[0], n2);
//  for (int p = 0; p < nsites; p++)
//  for (int q = 0; q < nsites; q++)
//    f[p*nsites+q] += one_int[p*nsites+q];
//  double energy = 0.5 * pcps::ddot(nsites*nsites, &f[0], 1, &rdm1[0], 1);
//  for (int p = 0; p < nsites; p++)
//  for (int q = 0; q < nsites; q++)
//    energy += 0.5 * rdm1[p*nsites+q] * one_int[p*nsites+q];
//
//  // Print the energy
//  std::cout << boost::format("free particle energy = %20.12f") % energy << std::endl;
//  std::cout << std::endl;
//
//  // Print the solution's 1-body RDM in the AO basis
//  std::cout << boost::format("1-body RDM:") << std::endl;
//  std::cout << std::endl;
//  for (int r = 0; r < nsites; r++) {
//    for (int s = 0; s < nsites; s++)
//      std::cout << boost::format("%12.5f") % rdm1[r*nsites+s];
//    std::cout << std::endl;
//  }
//  std::cout << std::endl;
//
//  // Print the 1-body RDM's trace
//  double trace = 0.0;
//  for (int r = 0; r < nsites; r++)
//    trace += rdm1[r*nsites+r];
//  std::cout << boost::format("trace = %.5f") % trace << std::endl;
//  std::cout << std::endl;
//
//  // Print the solution's MO coefficients
//  std::cout << boost::format("MO coefficients:") << std::endl;
//  std::cout << std::endl;
//  for (int i = 0; i < nsites; i++) {
//    if (i < nparticles)
//      std::cout << boost::format("%4s  ") % "occ";
//    else
//      std::cout << boost::format("%4s  ") % "vir";
//    for (int s = 0; s < nsites; s++)
//      std::cout << boost::format("%12.5f") % c[i*nsites+s];
//    std::cout << std::endl;
//  }
//  std::cout << std::endl;
//
//  // save 1-body RDM to disk
//  {
//    std::string file_path = userinp.workdir() + "rdm1.dat";
//    std::ofstream outfile(file_path.c_str(), std::ios::out | std::ios::binary);
//    if (!outfile.is_open())
//      throw pcps::Exception( (boost::format("failed to open file \"%s\" for saving 1-body RDM") % file_path).str() );
//    outfile.write( (const char *)(&rdm1[0]), rdm1.size() * sizeof(double) );
//    outfile.close();
//  }
//  
//  // Print the elapsed time
//  std::cout << boost::format("time to compute 1-body RDM for the free particle wavefunction: %.4f sec") % sw.elapsed_seconds() << std::endl;
//  std::cout << std::endl;
//}
